Spontaneous dimerization in the spin-1 bilinear-biquadratic Heisenberg model on a 

honeycomb lattice 
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Within the hnear flavor- wave theory, we show that, due to the quantum order- by-disorder mech- 
anism, the spin-1 bilinear-biquadratic Heisenberg model on a honeycomb lattice can spontaneously 
develop a columnar dimer order with a non-bipartite structure. The low-lying excitations above this 
novel ground state form several flat bands separated by nonzero energy gaps. Our results suggest 
that the quantum phase transition separating this dimerized phase with the nearby Neel-order phase 
may be of first order. 
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Quantum spin systems on various kinds of lattices have 
provided a wide playground for the search of novel quan- 
tum states and quantum phase transitions. It has been 
proposed that, in some spin- 1/2 systems on a square lat- 
tice, quantum paramagnetic phases with valence-bond- 
solid order can emerge and their transitions to the anti- 
ferromagnetic Neel states may belong to the dcconfincd 
quantum criticality.^'^ On a honeycomb lattice, even the 
more fascinating quantum spin-liquid state has recently 
been proposed and even observed in the half-filled Hub- 
bard model and related frustrated spin-1/2 Heisenberg 
mo dels. 

Rich physics is also expected in quantum spin systems 
with larger spin moments. Prominent examples include 
the spin-1 bilinear-biquadratic (BLBQ) model described 
by the Hamiltonian 



H 



E 



(1) 



with a bilinear exchange coupling J = cos 9 and a bi- 
quadratic exchange coupling K = sinO, where the param- 
eter 9 controls the ratio of these two couplings. It has 
been shown that quantum phases with either collinear 
or noncollinear nematic order can appear in this spin 
system. ^'^ In these nematic phases, dipolar spin order 
parameters vanish, (S) = 0, but spin rotation symme- 
try is spontaneously broken due to nonzero quadrupo- 
lar spin expectation values, {S" S/^ + S^^ S") ^ (a, 
f3 = X, y, z). Here we focus on the parameter region with 
tt/4: < 6 < tt/2 (i.e., K > J > 0), where the noncollinear 
nematic states described by mutually orthogonal nematic 
directors are expected at the mean-field level. When 
the BLBQ model is defined on a triangular lattice, in 
which each site has 6 neighbors, the ground state is found 
to posses a three-sublattice nematic order. ^'^ The ne- 
matic directors on the three sublattices A, B, and C of 
the triangular lattice arc orthogonal to each other (say, 
along the x, y, and z directions, respectively). On the 
other hand, if this model is placed on a square lattice, be- 
cause of the smaller coordination number for this lattice 
structure, there are not enough constraints to uniquely 
determine a set of mutually orthogonal nematic directors. 



Therefore, the expected noncollinear nematic state will 
form a highly degenerate ground-state manifold. Usu- 
ally, such macroscopic degeneracy at the variational level 
can be lifted by quantum effects, and a unique ground 
state will be selected by the order-by-disorder mecha- 
nism. Following this reasoning, the ground state for the 
BLBQ model with J = K on a square lattice has been 
analyzed recently by means of a semiclassical fiavor-wave 
theory and exact diagonalizations.^'' It is found that in- 
stead of the naive two-sublattice state, the ground state 
develops an unexpected three-sublattice long-range or- 
der, even though this system is defined on a bipartite 
lattice and with only nearest-neighbor interactions. 

Since the honeycomb lattice has an even smaller coor- 
dination number, one may wonder whether it can support 
ground states quite different from those on a square lat- 
tice. A possible candidate of the quantum phase for the 
BLBQ model on a honeycomb lattice with tt/A < 9 < Tr/2 
(i.e., K > J > 0) has been proposed. By employing 
the tensor renormalization group method^^'^'^ under the 
assumption that the ground-state wave function can be 
described by a periodic tensor network with elementary 
hexagon as the unit cell, the ground state is found to have 
plaquette valance-bond-solid (PVBS) order. This PVBS 
state breaks the lattice translational symmetry but pre- 
serves the spin SU(2) symmetry. However, due to the 
assumed periodicity for their variational states, the pos- 
sibilities for the ground states to display more compli- 
cated structures arc missed in their exploration. Hence 
their results may not be conclusive. 

In the present work, we determine the nature of the 
ground state of the BLBQ model in Eq. (1) with K > 
J > on a honeycomb lattice using the linear flavor- 
wave (LFW) theory.^''''^''^'-'^^ Previous studies on the cases 
of the square lattice have proved the validity and success 
of this approach. In Rcf. 10, the conclusion of the 
LFW analysis for the square-lattice case was shown to 
be supported by the exact diagonalization calculations 
for systems of finite sizes. Besides, the results within 
the LFW theory for the SU(4) Heisenberg model pro- 
vide useful guides for numerical methods to uncover the 
intricate dimerized structure and color ordering of the 
ground state. ^® For the present case of the spin-1 BLBQ 



model on a honeycomb lattice, following the reasoning of 
these works, we find that the ground states exhibit exotic 
columnar dimer order [see Fig. 1(c)] and break sponta- 
neously both the spin SU(2) symmetry and the lattice 
translation symmetry. Our results arise from the order- 
by-disorder mechanism, where the ground states are se- 
lected among the degenerate manifold of the noncollincar 
ncmatic states by minimizing the zero-point energies of 
quantum fluctuations. We note that the unit cell of the 
resulting ground-state pattern is quite large (consisting 
of 18 lattice sites), while the underlying lattice is bipar- 
tite. Within our LFW analysis, the flavor-wave excita- 
tions in the present systems are found to be localized, 
and the low-lying modes form several flat bands sepa- 
rated by nonzero energy gaps. This observation is rather 
different from the usual cases with spontaneous SU(2) 
symmetry breaking. Moreover, it is found that the gap- 
ful ground state with columnar dimer order remains even 
in the limit oi J = K (or 9 = 7r/4), at which a direct tran- 
sition to the nearby antiferromagnetic phase for J > K 
is anticipated.^'^ This suggests that the quantum phase 
transition aX J = K, which separates these two phases 
with distinct types of long-range order, may be of first or- 
der. The implication of our results on generalized models 
is discussed at the end of this Rapid Communication. 

The LFW theory starts from representing the model 
in Eq. (1) in terms of three-flavor Schwinger bosons ai^a 
under the local constraint V at q, = 1, 



H 



{K - J)AIA,, + [K ^ J) . (2) 



Here we define two bond operators, Xij = X^a ' 



J2a ' 



The Schwinger bosons 



and Ay- 

(with a =: X, y, z) create three time-reversal-invariant 
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-1)), 

-1)), and \z) = \s, = 0). 
In terms of these bosons, the spin operators become 



local basis states, \x) 

\y) = 71 (1^^ = 1) + 



''^p f^ap-yoJi ptti^^. The first step is the mean- 



field analysis based on a site-factorized variational wave 
function. At this mean-field level, the Schwinger-boson 



operators = (a^.^, ai 



;) are replaced by a (com- 



plex) three-component vector d^, and the energy of the 
nearest-neighbor bond (j, j) is minimal when the two vec- 
tors di and dj are mutually orthogonal. ^'^ The mean-field 
ground state configuration is highly degenerate on a hon- 
eycomb lattice, and as can be seen from Eq. (2), the asso- 
ciated mean-field energy per bond is [K—J). This macro- 
scopic degeneracy can be lifted when quantum fiuctua- 
tions above each mean-field state are included. Within 
the LFW analysis, the leading quantum corrections to 
the mean-field energy of the considered variational state 
come from the zero-point energy of the LFW Hamilto- 
nian. Therefore, the configuration with the lowest zero- 
point energy will be picked out as the true ground state. 
This is in essence a quantum order-by-disorder selection 
mechanism. 
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FIG. 1: (Color online) Schematic representation of three pos- 
sible candidates for the ground state: (a) two-sublattice state 
and (b) staggered dimer state, both of which have higher LFW 
energies; (c) columnar dimer state selected by quantum fluc- 
tuations within the LFW theory. Here the three mutually 
orthogonal vectors di in the mean-field analysis are denoted 
by the unit vectors along the x, y, and z directions, and are 
associated with different colors. Strong bonds with the lowest 
zero-point energy are denoted by thick lines. The remaining 
weak bonds are depicted by thin lines. 



The derivation of the LFW Hamiltonian proceeds as 
follows. For a given configuration of the variational state, 
when the two classical vectors on a nearest-neighbor bond 
(i, j) are, say, = x and dj = y (named as the XY 
bond in the following), we approximate semiclassically 
their Schwinger-boson operators by ~ (1, Ui^y, Ui^z) 
and SLj ~ io,j.x, 1, o,j.z)- The operators ai^y and Ui^z 
on site i {aj^x and aj^z on site j) describe the quantum 
fluctuations around the classical vectors, and they play 
the role of the Holstcin-Primakoff bosons in the usual 
spin- wave theory. Therefore, up to the linear order in 
these operators, the bond operators become 



Xij 



A,; 



(3) 



Substituting them into the Hamiltonian in Eq. (2), we 
obtain the desired quadratic LFW Hamiltonian. 

We note that the LFW Hamiltonian in general consists 
of a sum of independent parts that describe the motion 
of bosons on certain connected clusters. For example, for 
a given nearest-neighbor bond (say, the XY bond), when 
all of the d vectors of its surrounding sites are orthogonal 
(say, d; = z) to both d vectors on that bond, the 2 bosons 
(say. the x and y components of the Schwinger bosons) 
within that bond cannot move to neighboring sites. That 
is, the motion of these bosons becomes decoupled from 
the surrounding of that bond, and it can be described by 
a 2-site Hamiltonian. In the following, such 2-site clus- 
ters are dubbed as strong bonds and denoted by thick 
lines in Figs. 1(b) and 1(c). On the other hand, the re- 
maining bonds depicted by thin lines can be linked to 
from larger cluster and they are termed as weak bonds. 
As discussed below, we find that the zero-point energy 
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(and therefore the ground-state energy) is minimized for 
those strong bonds. Therefore, we expect that the states 
with the lowest zero-point energy should be the ones that 
contain as many strong bonds as possible. Such a "max- 
imum strong bond rule," which has been noted in other 
context,^® forms our main guiding principle in searching 
for possible candidates of the ground state. Based on 
this observation, in addition to the two-sublatticc state 
[Fig. 1(a)], which is the naive ground-state configuration 
on the present bipartite lattice, we consider two more 
configurations. They are the staggered dimer state and 
the columnar dimer state [Figs. 1(b) and Figs. 1(c), re- 
spectively] , both of which contain the maximum number 
of strong bonds per elementary hexagon. Nevertheless, 
as discussed in the paragraph below Eq. (6), the zero- 
point energy for weak bonds in the columnar dimer state 
is lower than that in the staggered dimer state. Thus 
the columnar dimer state in Fig. 1(c) with a more com- 
plicated structure is selected by the zero-point quantum 
fluctuations within the LFW theory. 

Now we begin to derive the explicit expressions of the 
zero-point energies for the three configurations shown in 
Fig. 1. For the case of the two-sublattice state formed 
by XY bonds only, the z components of the Schwinger 
bosons play no role and the LFW Hamiltonian reduces 
to an effective model for the two-component boson mix- 
ture on a single connected cluster of the whole honey- 
comb lattice. By means of the Bogoliubov transforma- 
tion, the excitation spectrum of the resulting quadratic 
bosonic Hamiltonian can be easily obtained. This re- 
sults in the following ground-state energy per site (i.e., 
the sum of the mean-field energy density and the zero- 
point contribution), Eq = -^^^i^ik + ^2k.)- Here 
N is the number of lattice sites, k runs over the first 
Brillouin zone of the honeycomb lattice, and Ai.2(k) = 

£0 ± |A2k|)' - jAikj' with eo = ^K, Aik = f J7k, 
and A2k = ^{K — J)"fu.. The coordination number z = 3 
for the honeycomb lattice, and 7k = 7 ^^^'^ , where S 
runs over the three nearest-neighbor vectors, 61 = (1,0), 



ground-state energy can be found as 



'52 = (-i^), and<53 = (-i,-^). 

Unlike the two-sublattice state, both the staggered 
dimer state and the columnar dimer state have strong 
bonds. Besides, their weak bonds form many one- 
dimensional zig-zag chains or 6-site clusters (hexagonal 
loops) as shown in Fig. 1(b) and Fig. 1(c), respectively. 
Since the bosons on these weak bonds are decoupled from 
those on the strong bonds, these bosons can be described 
by an effective Hamiltonian on a one-dimensional chain 
of 2Nc bonds, where Nc 00 and Nc ^ 3 for the stag- 
gered dimer state and the columnar dimer state, respec- 
tively. That is, the fiavor-wave excitations are well local- 
ized within those chains or clusters. Because the LFW 
Hamiltonians on clusters of different sizes give distinct 
contributions to the zero-point energy, they need to be 
considered separately. For a given strong bond, by diag- 
onalizing the corresponding 2-site LFW Hamiltonian, its 



s-bon 



d = ^K{K - J) 



(4) 



Similarly, for excitations localized within a closed chain 
of 2Nc bonds, diagonalization of the corresponding one- 
dimensional Hamiltonian results in the following expres- 
sion for the ground-state energy per weak bond. 



-bond 



(iVc 



(5) 



Here the onc-dimcnsional momentum k =^ with n 



l,...,iVe, Ai,2(fc) = \/ £0±|A 



|Aife|2 with eo 



K, jAifcl = J\ cos(|)|, and |A2fe| = {K - J)\ cos(|)|. As 
seen from Figs. 1(b) and 1(c), in both the staggered dimer 
state and the columnar dimer state, the three nearest- 
neighbor bonds for each lattice site contain one strong 
and two weak bonds. Therefore, from Eqs. (4) and (5), 
the expression of the ground-state energy per site for both 
cases becomes 



1 



(6) 



Thus the difference in energies between the staggered 
dimer state and the columnar dimer state comes from 
the contribution of the weak bonds. From Eq. (5), one 
can show that the energy per weak bond is an increasing 
function of Nc- In other words, the localized fiavor-wave 
excitations residing on shorter chains will give a smaller 
contribution in Eq. (6). As mentioned above, Nc ~ 3 
for the columnar dimer state and — >■ c» for the stag- 
gered dimer state. Hence, we conclude that the former 
has a lower energy. Note that on a honeycomb lattice, 
the shortest closed loop of weak bonds is nothing but the 
hexagonal one with Nc = 3. Thus we conclude that the 
columnar dimer state should be the true ground state 
among the degenerate manifold. 

The ground-state energy densities for the three types of 
states are presented in Fig. 2. We find that their energies 
become degenerate and approach to 3/2 in the — > 7r/2 
limit. This comes from the fact that the mean-field state 
becomes an exact eigenstate at = 7r/2 (or J = and 
K = 1) and quantum fluctuations play no role here. Thus 
the mean- field energy per site Eq = ^{K — J) = 3/2 
is nothing but the exact energy eigenvalue. Moreover, 
our Fig. 2 shows that the columnar dimer state does 
have the lowest energy in the whole parameter regime 
7r/4 < 6* = tan-i(K/J) < 7r/2. We note that the PVBS 
states proposed in Ref. 11 lie outside the degenerate man- 
ifold of the noncollincar ncmatic states and thus are not 
considered within the employed approach. However, we 
can compare their variational energies with those dis- 
cussed in the present work. It is found that the ground- 
state energies of our columnar dimer states are also lower 
than those of the PVBS states. For instance, at the 
SU(3) point with e = 7r/4 (or J = K = 1/a/2), our 
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FIG. 2: (Color online) Ground state energy densities (per 
site) for the three states shown in Fig. 1 as a function of 
6 — twa~^ [K / J) . The dash-dotted, dashed, and solid lines are 
for the two-sublattice, the staggered dimer, and the columnar 
dimer states, respectively. Inset: Energy difference (per site) 

ABg = Sg, staggered dimer — Sg, column dimer between the Stag- 
gered dimer and the columnar dimer states as a function of 

e. 



values of Eq for the two-sublattice state, the staggered 
dimer state, and the eolumnar dimer state are 0.8381, 
V2/7r (~ 0.4502), and (~ 0.4082), respectively. 

However, the energy density at = 7r/4 for the PVBS 
state is about 0.536, which is higher than the values for 
both of the dimerized states. Thus our calculations show 
that the peculiar dimerized ground states that sponta- 
neously break both spatial and spin SU(2) symmetries 
can in fact be energetically more favorable. We should 
stress that since the energies obtained within the LFW 
theory are not variational, the above comparison in en- 
ergy may not be used to exclude the PVBS states as the 
true ground state. Nevertheless, the present investiga- 
tion indicates that our dimerized states should at least 
be possible candidates. 

To provide guides for future numerical investigations, 
some comments are discussed below. People usually per- 
form exact diagonalization calculations for systems of fi- 
nite sizes to obtain unbiased results and then provide 
numerical verification for analytical proposals. Since the 
unit cell of our eolumnar dimer state contains 18 lattice 
sites, large systems with sizes of several unit cells should 
be considered in order to deliver meaningful results com- 
plementary to our present work. But the required com- 
putational resources may make such calculations unlikely. 
In this regard, numerical variational approaches, such as 
the variational Monte Carlo method and the tensor net- 
work method, may be more promising, because their size 
limitation is less severe. Even in this case, our results 
show that one has to use the variational states compati- 
ble with the spatial structure of our columnar dimer state 
to determine convincingly the true ground state. 

In addition to determining the nature of the ground 



state within parameter regime 7r/4 < 9 < tt/2, our re- 
sults also help to characterize the types of phase transi- 
tions out of this phase. It is known that at the mean- 
field level, a phase transition to the nearby ferromagnetic 
phase occurs at 6* = 7r/2.^'^ At this value of 6, we find 
that the energy of the columnar dimer state {Eq —3/2) 
is identical to that of the ferromagnetic state [see Eq. (1) 
with Si = 1] . Because these two states belong to different 
subspaces of the total spin, this agreement in energy in- 
dicates a level crossing and thus a first-order transition at 
6 = tt/2. On the other hand, the mean-field theory pre- 
dicts a direct transition at 9 = tt/A from the noneoUinear 
nematic to the antiferromagnetic phases. ^'^ In our work, 
the effects of quantum fluctuations are taken into account 
under the LFW analysis. We find that there exist only 
the localized fiavor-wave excitations above the columnar 
dimer state and these excitations form fiat bands sep- 
arated by nonzero energy gaps. Moreover, the gapful 
ground state with nonvanishing eolumnar dimer order 
persists even in the limit of 6* = 7r/4 (or J = K). This 
suggests that the quantum phase transition at J = A', 
which separates these two phases with distinct types of 
long-range order, should be of first order. 

Finally, the implication of our results on generalized 
models is discussed. Unlike its counterparts on other 
lattices with larger coordination numbers z [i.e., the 
three-sublattice antiferronematic state on a triangular 
lattice^'^ [z = 6) and the three-sublattice stripelike state 
on a square lattice^*^ (z = 4)], the columnar dimer state 
on a honeycomb lattice (z = 3) has a solidlike pat- 
tern and supports only localized fiavor-wave excitations 
within the LFW analysis. In addition, a conclusion 
similar to ours is reached for an SU(4)-symmetric spin 
model with 4 states in the local Hilbert space. From 
all these findings, one may conclude that for generalized 
models with more local states [say, the SO(ri) bilinear- 
biquadratic model^®"^°] and/or on two-dimensional lat- 
tices with smaller z, ground states with solid-like struc- 
tures will be stabilized due to the same order-by-disorder 
mechanism. While this observation is based on the LFW 
analysis, we believe that the qualitative picture will not 
be modified even if higher quantum corrections are in- 
cluded. 

To summarize, by employing the LFW theory, we show 
that the columnar dimer state is the ground state of 
the spin-1 BLBQ model on a honeycomb lattice with 
7r/4 < 9 < tt/2. We stress that our investigations are 
not of purely academic interest. In fact, aX = n/A (i.e., 
K — J > 0), the present model possesses an enlarged 
SU(3) symmetry and can be considered as an effective 
model for the Mott-insulating state of three-fiavor cold 
fermions with one particle per lattice site.^° Therefore, 
our conclusions at this SU(3) point may be examined 
experimentally for such cold fermions on a honeycomb 
optical lattice. 
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